pacman::p_load(tidyverse, data.table, broom, openxlsx, kableExtra)
rm(list=ls())
################################################################################

######################## PARAMETERS 

###Spatial units
df_u = setDT(read.xlsx('aux_params/parameters.xlsx', sheet='spatial_units'))
df_u1 = setDT(read.xlsx('aux_params/parameters.xlsx', sheet='spatial_units'))[,commodity:='beef']
df_u2 = setDT(read.xlsx('aux_params/parameters.xlsx', sheet='spatial_units'))[,commodity:='crops']
df_u_c = rbind(df_u1,df_u2)[,c('spatial_id_name','commodity')]

###Emissions
params_eLUC = setDT(read.xlsx('aux_params/parameters.xlsx',  sheet = "emissions_luc"))[,list(spatial_id_name,tCO2_per_ha_flow)]
colnames(params_eLUC) = c('spatial_id_name','eLUC')
params_eNLUC = setDT(read.xlsx('aux_params/parameters.xlsx',  sheet = "emissions_nluc"))[,list(spatial_id_name,beef_nluc_per_t,crops_nluc_per_t)]
colnames(params_eNLUC) = c('spatial_id_name','beef','crops')
params_eQ = setDT(read.xlsx('aux_params/parameters.xlsx',  sheet = "emissions_total"))[,list(spatial_id_name,beef_total_flow_per_t,crops_total_flow_per_t)]
colnames(params_eQ) = c('spatial_id_name','beef','crops')

###SCC
df = read.xlsx('aux_params/scc.xlsx', sheet="baseline_ic")
SCC = max(df$SCC)

###Demand
df = setDT(read.xlsx('aux_params/parameters.xlsx',  sheet = "indicator_in"))
df$origin_nation = 'n'
df[source_a==1,]$origin_nation = 'ARGENTINA'
df[source_b==1,]$origin_nation = 'BRAZIL'
df = df[,list(spatial_id_name,origin_nation)]
params_in = df

df = setDT(read.xlsx('aux_params/parameters.xlsx',  sheet = "demand_j"))
n <- df$variable
df <- as.data.frame(t(df[,-1]))
colnames(df) <- n
df$destination <- factor(row.names(df))
params_demand_j = df

df = setDT(read.xlsx('aux_params/parameters.xlsx',  sheet = "a_j_c"))
setnames(df, old=c('product'), new=c('commodity'))
df <- melt(df, id.vars = c("commodity"), variable.name = "destination", value.name = "a_j_c")
params_demand_jc = df

df = setDT(read.xlsx('aux_params/parameters.xlsx',  sheet = "a_nj_c"))
setnames(df, old=c('origin_nation','product'), new=c('origin_nation','commodity'))
df <- melt(df, id.vars = c("origin_nation","commodity"), variable.name = "destination", value.name = "a_nj_c")
params_demand_njc = df

df = setDT(read.xlsx('aux_params/parameters.xlsx',  sheet = "a_ij_c"))
setnames(df, old=c('spatial_id_name','product'), new=c('spatial_id_name','commodity'))
df <- melt(df, id.vars = c("spatial_id_name","commodity"), variable.name = "destination", value.name = "a_ij_c")
params_demand_ijc = df

###MP
df = setDT(read.xlsx('aux_params/parameters.xlsx',  sheet = "MP"))
setnames(df, old=c('mp_beef','mp_crops'), new=c('beef','crops'))
params_mp <- melt(df, id.vars = c("spatial_id_name"), variable.name = "commodity", value.name = "mp")



######################## WELFARE 

base_spec = 'lf'
cf_list = c('tariff','dt','ur','certr','certc')
cf_list_fb = append(cf_list,'fb')

cf_spec_n=1
for(cf_spec in cf_list_fb){
  
  #################### MARKET OUTCOMES
  
  #Production
  df = setDT(read.csv(file = paste0('results_ic/',base_spec,'/Q.csv'), header=FALSE, col.names = c('beef','crops')))
  df$spatial_id_name = df_u$spatial_id_name
  df_0 <- melt(df, id.vars = c("spatial_id_name"), variable.name = "commodity", value.name = "level0")
  
  df =  setDT(read.csv(file = paste0('results_ic/',cf_spec,'/Q.csv'), header=FALSE, col.names = c('beef','crops')))
  df$spatial_id_name = df_u$spatial_id_name
  df_1 <- melt(df, id.vars = c("spatial_id_name"), variable.name = "commodity", value.name = "level1")
  
  df_Q = merge(df_0,df_1, by = c("spatial_id_name","commodity"))[, delta:=level1-level0]
  setnames(df_Q, old=c('level0','level1','delta'), new=c('Q0','Q1','deltaQ'))
  
  #Consumption
  df = setDT(read.csv(file = paste0('results_ic/',base_spec,'/C_ij.csv'), header=FALSE, col.names = c('SA','EU','AS','ROW')))
  df$spatial_id_name = df_u_c$spatial_id_name
  df$commodity = df_u_c$commodity
  df_0 <- melt(df, id.vars = c("commodity","spatial_id_name"), variable.name = "destination", value.name = "level0")
  
  df = setDT(read.csv(file = paste0('results_ic/',cf_spec,'/C_ij.csv'), header=FALSE, col.names = c('SA','EU','AS','ROW')))
  df$spatial_id_name = df_u_c$spatial_id_name
  df$commodity = df_u_c$commodity
  df_1 <- melt(df, id.vars = c("commodity","spatial_id_name"), variable.name = "destination", value.name = "level1")
  
  df_C = merge(df_0,df_1, by = c("spatial_id_name","destination","commodity"))[, delta:=level1-level0]
  setnames(df_C, old=c('level0','level1','delta'), new=c('C0','C1','deltaC'))
  
  #Land
  df = setDT(read.csv(file = paste0('results_ic/',base_spec,'/L.csv'), header=FALSE, col.names = c('beef','crops','forest')))
  df$spatial_id_name = df_u$spatial_id_name
  df_0 <- melt(df, id.vars = c("spatial_id_name"), variable.name = "commodity", value.name = "level0")
  
  df = setDT(read.csv(file = paste0('results_ic/',cf_spec,'/L.csv'), header=FALSE, col.names = c('beef','crops','forest')))
  df$spatial_id_name = df_u$spatial_id_name
  df_1 <- melt(df, id.vars = c("spatial_id_name"), variable.name = "commodity", value.name = "level1")
  
  df_L = merge(df_0,df_1, by = c("spatial_id_name","commodity"))[, delta:=level1-level0]
  setnames(df_L, old=c('level0','level1','delta'), new=c('L0','L1','deltaL'))
  
  #Prices - farm-gate - pre-tax and after-tax
  df = setDT(read.csv(file = paste0('results_ic/',base_spec,'/p_i_c_pretax.csv'), header=FALSE, col.names = c('beef','crops')))
  df$spatial_id_name = df_u$spatial_id_name
  dfx <- melt(df, id.vars = c("spatial_id_name"), variable.name = "commodity", value.name = "level0")
  df = setDT(read.csv(file = paste0('results_ic/',base_spec,'/p_i_c.csv'), header=FALSE, col.names = c('beef','crops')))
  df$spatial_id_name = df_u$spatial_id_name
  dfy <- melt(df, id.vars = c("spatial_id_name"), variable.name = "commodity", value.name = "level0_at")
  df_p_i_0 = merge(dfx,dfy, by= c("spatial_id_name","commodity"))
  
  df = setDT(read.csv(file = paste0('results_ic/',cf_spec,'/p_i_c_pretax.csv'), header=FALSE, col.names = c('beef','crops')))
  df$spatial_id_name = df_u$spatial_id_name
  dfx <- melt(df, id.vars = c("spatial_id_name"), variable.name = "commodity", value.name = "level1")
  df = setDT(read.csv(file = paste0('results_ic/',cf_spec,'/p_i_c.csv'), header=FALSE, col.names = c('beef','crops')))
  df$spatial_id_name = df_u$spatial_id_name
  dfy <- melt(df, id.vars = c("spatial_id_name"), variable.name = "commodity", value.name = "level1_at")
  df_p_i_1 = merge(dfx,dfy, by= c("spatial_id_name","commodity"))
  
  df_p_i = merge(df_p_i_0,df_p_i_1, by= c("spatial_id_name","commodity"))[, c('delta','delta_at'):=list(level1-level0,level1_at-level0_at)]
  setnames(df_p_i, old=c('level0','level1','delta','level0_at','level1_at','delta_at'), new=c('p_i_0','p_i_1','deltap_i','p_at_i_0','p_at_i_1','deltap_at_i'))
  
  #Prices - consumer - pre-tax and after-tax
  df = setDT(read.csv(file = paste0('results_ic/',base_spec,'/p_ij_c_pretax.csv'), header=FALSE, col.names = c('SA','EU','AS','ROW')))
  df$spatial_id_name = df_u_c$spatial_id_name
  df$commodity = df_u_c$commodity
  dfx <- melt(df, id.vars = c("commodity","spatial_id_name"), variable.name = "destination", value.name = "level0")
  df = setDT(read.csv(file = paste0('results_ic/',base_spec,'/p_ij_c.csv'), header=FALSE, col.names = c('SA','EU','AS','ROW')))
  df$spatial_id_name = df_u_c$spatial_id_name
  df$commodity = df_u_c$commodity
  dfy <- melt(df, id.vars = c("commodity","spatial_id_name"), variable.name = "destination", value.name = "level0_at")
  df_p_ij_0 = merge(dfx,dfy, c("spatial_id_name","destination","commodity"))
  
  df = setDT(read.csv(file = paste0('results_ic/',cf_spec,'/p_ij_c_pretax.csv'), header=FALSE, col.names = c('SA','EU','AS','ROW')))
  df$spatial_id_name = df_u_c$spatial_id_name
  df$commodity = df_u_c$commodity
  dfx <- melt(df, id.vars = c("commodity","spatial_id_name"), variable.name = "destination", value.name = "level1")
  df = setDT(read.csv(file = paste0('results_ic/',cf_spec,'/p_ij_c.csv'), header=FALSE, col.names = c('SA','EU','AS','ROW')))
  df$spatial_id_name = df_u_c$spatial_id_name
  df$commodity = df_u_c$commodity
  dfy <- melt(df, id.vars = c("commodity","spatial_id_name"), variable.name = "destination", value.name = "level1_at")
  df_p_ij_1 = merge(dfx,dfy, c("spatial_id_name","destination","commodity"))
  
  df_p_ij = merge(df_p_ij_0, df_p_ij_1, by=c("spatial_id_name","destination","commodity"))[, c('delta','delta_at'):=list(level1-level0,level1_at-level0_at)]
  setnames(df_p_ij, old=c('level0','level1','delta','level0_at','level1_at','delta_at'), new=c('p_ij_0','p_ij_1','deltap_ij','p_at_ij_0','p_at_ij_1','deltap_at_ij'))
  
  #Markdowns
  df = setDT(read.csv(file = paste0('results_ic/',base_spec,'/mu_i_c.csv'), header=FALSE, col.names = c('beef','crops')))
  df$spatial_id_name = df_u$spatial_id_name
  df_0 <- melt(df, id.vars = c("spatial_id_name"), variable.name = "commodity", value.name = "level0")
  
  df =  setDT(read.csv(file = paste0('results_ic/',cf_spec,'/mu_i_c.csv'), header=FALSE, col.names = c('beef','crops')))
  df$spatial_id_name = df_u$spatial_id_name
  df_1 <- melt(df, id.vars = c("spatial_id_name"), variable.name = "commodity", value.name = "level1")
  
  df_mu = merge(df_0,df_1, by = c("spatial_id_name","commodity"))[, delta:=level1-level0]
  setnames(df_mu, old=c('level0','level1','delta'), new=c('mu_0','mu_1','deltamu'))

  #Intermediary prices
  df=merge(df_p_i,df_mu, by=c('spatial_id_name','commodity'))
  df=merge(df,params_mp, by=c('spatial_id_name','commodity'))
  df=df[, c('p_ibar_0','p_ibar_1'):=list( (p_i_0/mu_0)/mp , (p_i_1/mu_1)/mp ) ][,deltap_i_bar:=p_ibar_1-p_ibar_0]
  df_p_ibar=df[, list(spatial_id_name, commodity, mu_0, mu_1, deltamu, p_ibar_0, p_ibar_1,deltap_i_bar)]
  
  
  
  #################### EMISSIONS
  
  #LUC emissions
  df_forest = df_L[commodity=='forest',]
  df_total  = df_L[,list(spatial_id_name,L0,L1)][, lapply(.SD, sum, na.rm=TRUE), by=list(spatial_id_name)][,L0_total:=L0][,L1_total:=L1][,-c('L0','L1')]
  df = merge(df_forest,df_total,by=c('spatial_id_name'))
  df = merge(params_eLUC,df, by=c('spatial_id_name'))[, c('E_L0','E_L1','deltaE_L'):=list(eLUC*(L0_total-L0),eLUC*(L1_total-L1),-eLUC*deltaL)]
  deltaE_L = sum(df$deltaE_L)
  E_L0 = sum(df$E_L0)

  #NLUC emissions
  df <- melt(params_eNLUC, id.vars = c("spatial_id_name"), variable.name = "commodity", value.name = "eNLUC")
  df = merge(df,df_Q, by=c('spatial_id_name','commodity'))[, c('E_NL0','E_NL1','deltaE_NL'):=list(eNLUC*Q0,eNLUC*Q1,eNLUC*deltaQ)]
  deltaE_NL = sum(df$deltaE_NL)
  E_NL0 = sum(df$E_NL0)

  #Total emissions
  deltaE = deltaE_L + deltaE_NL
  E0 = E_NL0 + E_L0
  deltaEpp = 100*deltaE/E0
  
  
  
  #################### GOVERNMENT SURPLUS
  
  #Downstream consumer taxes
  df = setDT(read.csv(file = paste0('results_ic/',base_spec,'/tax_ij.csv'), header=FALSE, col.names = c('SA','EU','AS','ROW')))
  df$spatial_id_name = df_u_c$spatial_id_name
  df$commodity = df_u_c$commodity
  df_tax_ij_0 <- melt(df, id.vars = c("commodity","spatial_id_name"), variable.name = "destination", value.name = "tax_ij_0")
  
  df = setDT(read.csv(file = paste0('results_ic/',cf_spec,'/tax_ij.csv'), header=FALSE, col.names = c('SA','EU','AS','ROW')))
  df$spatial_id_name = df_u_c$spatial_id_name
  df$commodity = df_u_c$commodity
  df_tax_ij_1 <- melt(df, id.vars = c("commodity","spatial_id_name"), variable.name = "destination", value.name = "tax_ij_1")
  
  df_tax_ij = merge(df_tax_ij_0,df_tax_ij_1, by = c("spatial_id_name","destination","commodity"))
  df = merge(df_C, df_tax_ij, by=c("commodity","spatial_id_name","destination"))[, c('GC1','GC0'):=list(C1*tax_ij_1, C0*tax_ij_0)][, deltaGC:=GC1-GC0]
  df = df[, list(commodity,destination,deltaGC)][, lapply(.SD, sum, na.rm=TRUE), by=list(commodity,destination)]
  if(cf_spec %in% c('dt','certr','certc','fb')){
    dfx = df
    dfy = df[, list(commodity,deltaGC)][, lapply(.SD, sum, na.rm=TRUE), by=list(commodity)]
    df = merge(dfx,dfy, by='commodity')
    df$deltaGC = (df$destination=='SA')*df$deltaGC.y
    df = df[, list(commodity,destination,deltaGC)] }
  deltaGC = sum(df$deltaGC)

  #Upstream output taxes
  df = setDT(read.csv(file = paste0('results_ic/',base_spec,'/tax_i.csv'), header=FALSE, col.names = c('beef','crops')))
  df$spatial_id_name = df_u$spatial_id_name
  df_tax_i_0 <- melt(setDT(df), id.vars = c("spatial_id_name"), variable.name = "commodity", value.name = "tax_i_0")
  
  df = setDT(read.csv(file = paste0('results_ic/',cf_spec,'/tax_i.csv'), header=FALSE, col.names = c('beef','crops')))
  df$spatial_id_name = df_u$spatial_id_name
  df_tax_i_1 <- melt(setDT(df), id.vars = c("spatial_id_name"), variable.name = "commodity", value.name = "tax_i_1")
  
  df_tax_i = merge(df_tax_i_0,df_tax_i_1, by = c("spatial_id_name","commodity"))
  df = merge(df_Q, df_tax_i, by=c("commodity","spatial_id_name"))[, c('GQ1','GQ0'):=list(Q1*tax_i_1, Q0*tax_i_0)][, deltaGQ:=GQ1-GQ0]
  df = df[, list(commodity,deltaGQ,GQ0,GQ1)][, lapply(.SD, sum, na.rm=TRUE), by=list(commodity)]
  deltaGQ = sum(df$deltaGQ)

  #Upstream land subsidies
  df = setDT(read.csv(file = paste0('results_ic/',base_spec,'/subsidy_i.csv'), header=FALSE, col.names = c('forest')))
  df$spatial_id_name = df_u$spatial_id_name
  df_subsidy_i_0 <- melt(df, id.vars = c("spatial_id_name"), variable.name = "commodity", value.name = "subsidy_i_0")
  
  df = setDT(read.csv(file = paste0('results_ic/',cf_spec,'/subsidy_i.csv'), header=FALSE, col.names = c('forest')))
  df$spatial_id_name = df_u$spatial_id_name
  df_subsidy_i_1 <- melt(df, id.vars = c("spatial_id_name"), variable.name = "commodity", value.name = "subsidy_i_1")
  
  df_subsidy_i = merge(df_subsidy_i_0,df_subsidy_i_1, by = c("spatial_id_name","commodity"))
  df = merge(df_L, df_subsidy_i, by=c("commodity","spatial_id_name"))[, c('GL1','GL0'):=list(-L1*subsidy_i_1,-L0*subsidy_i_0)][, deltaGL:=GL1-GL0]
  deltaGL = sum(df$deltaGL)
  GL0 =  sum(df$GL0)
  df_GL = df[, list(spatial_id_name,commodity,deltaGL,GL0,GL1)]
  
  #Total government surplus
  deltaG = deltaGC + deltaGQ + deltaGL
  
  
  
  #################### PRODUCER SURPLUS

  #Producer
  df = merge(df_Q,df_p_i, by=c("spatial_id_name","commodity"))[,c('PS0','PS1'):=list( Q0*(p_at_i_0), Q1*(p_at_i_1) )][, deltaPS:=PS1-PS0]
  df_PS_R = df[, list(commodity,PS0,deltaPS)][, lapply(.SD, sum, na.rm=TRUE), by=list(commodity)]
  df_PS_S = df_GL[, c('PS0','deltaPS'):=list(-GL0, -deltaGL)][, list(commodity,PS0,deltaPS)][, lapply(.SD, sum, na.rm=TRUE), by=list(commodity)]
  df = rbind(df_PS_R,df_PS_S)
  deltaPS = sum(df$deltaPS)
  PS0 = sum(df$PS0)
  deltaPSpp = 100*deltaPS/PS0

  #Intermediary
  df = merge(df_Q,df_p_ibar, by=c("spatial_id_name","commodity"))
  df = merge(df,params_mp, by=c("spatial_id_name","commodity"))
  df = merge(df,df_p_i, by=c("spatial_id_name","commodity"))
  df = df[, c('PSI0','PSI1'):=list( (p_ibar_0*mp-p_i_0)*Q0, (p_ibar_1*mp-p_i_1)*Q1 )][, deltaPSI:=PSI1-PSI0]
  deltaPSI = sum(df$deltaPSI)
  
  
  
  #################### CONSUMER SURPLUS
  
  eta_l = max(params_demand_j$eta_l)
  eta_m = max(params_demand_j$eta_m)
  eta_u = max(params_demand_j$eta_u)
  
  ####### Price indices 
  #Lower level
  df = merge(df_p_ij[,list(commodity, destination, spatial_id_name,p_at_ij_0,p_at_ij_1)], params_demand_ijc, by=c('commodity','destination','spatial_id_name'))
  df = merge(params_in, df, by='spatial_id_name')[, c('x0','x1'):=list(a_ij_c*(p_at_ij_0^(1-eta_l)),a_ij_c*(p_at_ij_1^(1-eta_l)))]
  df = df[, list(commodity,destination,origin_nation,x0,x1)][, lapply(.SD, sum, na.rm=TRUE), by=list(commodity, destination, origin_nation)]
  df = df[, c('P_nj_c0','P_nj_c1'):=list(x0^(1/(1-eta_l)),x1^(1/(1-eta_l)))][, delta:=P_nj_c1-P_nj_c0]
  #Middle level
  df = merge(df[,list(commodity, destination, origin_nation, P_nj_c0, P_nj_c1)], params_demand_njc, by=c('commodity','destination','origin_nation'))[, c('x0','x1'):=list(a_nj_c*(P_nj_c0^(1-eta_m)),a_nj_c*(P_nj_c1^(1-eta_m)))]
  df = df[, list(commodity,destination, x0, x1)][, lapply(.SD, sum, na.rm=TRUE), by=list(commodity, destination)]
  df = df[, c('P_j_c0','P_j_c1'):=list( x0^(1/(1-eta_m)),x1^(1/(1-eta_m)) )][, delta:=P_j_c1-P_j_c0]
  #Upper level
  df = merge(df[,list(commodity,destination, P_j_c0, P_j_c1)], params_demand_jc, by=c('commodity','destination'))[, c('x0','x1'):=list( a_j_c*(P_j_c0^(1-eta_u)), a_j_c*(P_j_c1^(1-eta_u)) ) ]
  df = df[, list(destination, x0, x1)][, lapply(.SD, sum, na.rm=TRUE), by=list(destination)]
  df = df[, c('P_j0','P_j1'):=list( x0^(1/(1-eta_u)),  x1^(1/(1-eta_u)) )]
  dfx = df[,list(destination,P_j0,P_j1)][, delta:=P_j1-P_j0]
  
  ####### Consumption aggregators 
  #Lower level
  df = merge(df_C[,list(commodity, destination, spatial_id_name,C0,C1)], params_demand_ijc, by=c('commodity','destination','spatial_id_name'))
  df = merge(params_in, df, by='spatial_id_name')[, c('x0','x1'):= list( (a_ij_c^(1/eta_l))*(C0^((eta_l-1)/eta_l)), (a_ij_c^(1/eta_l))*(C1^((eta_l-1)/eta_l)) ) ]
  df = df[, list(commodity,destination,origin_nation,x0,x1)][, lapply(.SD, sum, na.rm=TRUE), by=list(commodity, destination, origin_nation)]
  df = df[, c('C_nj_c0','C_nj_c1'):=list( x0^(eta_l/(eta_l-1)), x1^(eta_l/(eta_l-1)) )][, delta:=C_nj_c1-C_nj_c0]
  #Middle level
  df = merge(df[,list(commodity, destination, origin_nation, C_nj_c0, C_nj_c1)], params_demand_njc, by=c('commodity','destination','origin_nation'))[, c('x0','x1'):= list( (a_nj_c^(1/eta_m))*(C_nj_c0^((eta_m-1)/eta_m)), (a_nj_c^(1/eta_m))*(C_nj_c1^((eta_m-1)/eta_m)) )  ]
  df = df[, list(commodity,destination, x0, x1)][, lapply(.SD, sum, na.rm=TRUE), by=list(commodity, destination)]
  df = df[, c('C_j_c0','C_j_c1'):=list( x0^(eta_m/(eta_m-1)), x1^(eta_m/(eta_m-1)) )][, delta:=C_j_c1-C_j_c0]
  #Upper level
  df = merge(df[,list(commodity,destination, C_j_c0, C_j_c1)], params_demand_jc, by=c('commodity','destination'))[, c('x0','x1'):=list( (a_j_c^(1/eta_u))*(C_j_c0^((eta_u-1)/eta_u)), (a_j_c^(1/eta_u))*(C_j_c1^((eta_u-1)/eta_u)) )]
  df = df[, list(destination, x0, x1)][, lapply(.SD, sum, na.rm=TRUE), by=list(destination)]
  df = df[, c('U_j0','U_j1'):=list( x0^(eta_u/(eta_u-1)), x1^(eta_u/(eta_u-1)))]
  dfy = df[,list(destination,U_j0, U_j1)][, delta:=U_j1-U_j0]
  
  ####### Equivalent variation
  df = merge(dfx,dfy, by='destination')
  df_CS = merge(df, setDT(params_demand_j)[,list(destination, X_jt)], by='destination')[,x := P_j0*U_j1-X_jt]
  
  deltaCS = sum(df_CS$x)
  deltaCS_SA = df_CS[destination=='SA']$x
  deltaCS_nonSA = df_CS[destination=='EU']$x+df_CS[destination=='AS']$x+df_CS[destination=='ROW']$x

  CS0 = sum(params_demand_j$X_jt)
  CS0_SA =  df_CS[destination=='SA']$X_jt
  CS0_nonSA =  df_CS[destination=='EU']$X_jt+df_CS[destination=='AS']$X_jt+df_CS[destination=='ROW']$X_jt
  
  deltaCSpp = 100*deltaCS/CS0
  deltaCS_SApp = 100*deltaCS_SA/CS0_SA
  deltaCS_nonSApp = 100*deltaCS_nonSA/CS0_nonSA
  
  
  
  #################### WELFARE
  
  deltaW = deltaCS + deltaPS + deltaPSI + deltaG - SCC*deltaE
  tab = data.table(variable = c('W','CS','PS','E','CS_SA','CS_nonSA'),
                   c(deltaW, deltaCS, deltaPS, deltaE, deltaCS_SA, deltaCS_nonSA),
                   c(NA, deltaCSpp, deltaPSpp, deltaEpp, deltaCS_SApp, deltaCS_nonSApp))
  setnames(tab, c('variable','delta','deltapp'))
  n <- tab$variable
  tab <- as.data.frame(t(tab[,-1]))
  colnames(tab) <- n
  tab$variable <- factor(row.names(tab))
  tab$counterfactual = cf_spec
  
  if(cf_spec_n==1){tab_specs = tab}else{tab_specs = rbind(tab_specs,tab)}
  
  cf_spec_n=cf_spec_n+1
  
}


######################## SUMMARY TABLE

df = setDT(tab_specs)
dfx = cbind(df[variable=='delta', list(counterfactual, variable)],round(df[variable=='delta',     list(W,CS,PS,E,CS_SA,CS_nonSA)]/1000000,1))
dfy = cbind(df[variable=='deltapp', list(counterfactual, variable)],round(df[variable=='deltapp', list(W,CS,PS,E,CS_SA,CS_nonSA)],2))
df = rbind(dfx,dfy)

#Set row order - cf category
df$counterfactual_order = 1
df[counterfactual=='dt',]$counterfactual_order = 1
df[counterfactual=='certr',]$counterfactual_order = 2
df[counterfactual=='certc',]$counterfactual_order = 3
df[counterfactual=='tariff',]$counterfactual_order = 4
df[counterfactual=='ur',]$counterfactual_order = 5
df[counterfactual=='fb',]$counterfactual_order = 6
cf_labels = c('Untargeted',
              'Region-targeted',
              'County-targeted',
              'Unilateral EU tariff',
              'Forest subsidy')
cf_labels_spaced = c('Untargeted',' ',
                     'Region-targeted',' ',
                     'County-targeted',' ',
                     'Unilateral EU tariff',' ',
                     'Forest subsidy',' ')
df$variable_order = 1
df[variable=='delta',]$variable_order = 1
df[variable=='deltapp',]$variable_order = 2
df = df[order(variable_order,counterfactual_order)]

#Share of first-best efficiency gains
df_w = df[variable=='delta', list(counterfactual,W,E)]
df_w_fb = df_w[counterfactual=='fb', list(W,E)]
df_w = df_w[,list(W,E)]
df_w$W = round(100*df_w$W/df_w_fb$W,2)
df_w$E =  round(100*df_w$E/df_w_fb$E,2)
df_w_sharefb = df_w

#P.P. changes
df_w = df[variable=='deltapp',]
df_w$W_share = df_w_sharefb$W
df_w$E_share = df_w_sharefb$E
df_w_final=df_w[,list(E,PS,CS,CS_SA,CS_nonSA,E_share,W_share)]

#Final table
cf_label_fb = 'First-best'
cf_labels_final = append(cf_labels,cf_label_fb)
cf_labels_final_spaced = append(append(cf_labels_spaced,cf_label_fb),' ')
tab_aggregate = data.table(Counterfactual = cf_labels_final, df_w_final)
setnames(tab_aggregate, c('','$\\Delta$ E','$\\Delta$ PS','$\\Delta$ CS','$\\Delta$ CS$_{SA}$','$\\Delta$ CS$_{ROW}$','$\\Delta$ E','$\\Delta$ W'))
export_tab = knitr::kable(tab_aggregate, booktabs = T, align = "lccccccc", format = 'latex',linesep = c("","","","","",""),
                        caption=paste0('Decomposition of welfare impacts.'), label=paste0(' welfare - decomposition - final'), 
                        format.args = list(big.mark = ","), digits=2, escape=F) %>%
add_header_above(c(" " = 1, "equilibrium (percentage points)" = 5, "efficiency gains" = 2))  %>%
add_header_above(c(" " = 1, "Change relative to laissez-faire" = 5, "% of first-best" = 2))  %>%
kable_styling(latex_options = c("HOLD_position"))  %>% 
kableExtra::save_kable(export_tab, file = paste0("../../output/tables/table5_welfare.tex"))



